library(foreign)
library(psych)
library(Hmisc)
library(survey)
library(car)
require(nnet)
require(ggplot2)
require(reshape2)
library(psych)
library(lavaan)
library(stargazer)
library(nnet)
library(margins)
library(effects)


############This work is based on the 2016 wave of the European Social Survey
############Use the first line of the code to type the route to the ESS2016 file in your computer
############You can download the original file from https://www.europeansocialsurvey.org/download.html?file=ESS8e02_2&y=2016

############Select from the original dataset the 13 countries chosen and the individuals who voted in the past national elections 
dataset<-read.dta("C:\\ESS8e02.dta", convert.factors=F)
as.factor(dataset$cntry)
newdata <- dataset[ which(dataset$cntry =="AT" 
                          | dataset$cntry == "BE"
                          | dataset$cntry == "CH"
                          | dataset$cntry == "DE"
                          | dataset$cntry == "FI"
                          | dataset$cntry == "IT"
                          | dataset$cntry == "FR"
                          | dataset$cntry == "NL"
                          | dataset$cntry == "NO"
                          | dataset$cntry == "PL"
                          | dataset$cntry == "SE"
                          | dataset$cntry == "ES"
                          | dataset$cntry == "GB"), ]

newdata <- newdata[ which(newdata$vote==1), ]

########################Voting for populist and mainstream parties 
newdata$populist[newdata$prtvtbat== 3] <- 1
newdata$populist[newdata$prtvtcbe== 4|newdata$prtvtcbe== 7] <- 1
newdata$populist[newdata$prtvtfch== 1] <- 1
newdata$populist[newdata$prtvede1== 3] <- 1
newdata$populist[newdata$prtvede1== 6] <- 1
newdata$populist[newdata$prtvtdfi== 4] <- 1
newdata$populist[newdata$prtvtbit== 4|newdata$prtvtbit== 8|newdata$prtvtbit== 9] <- 1
newdata$populist[newdata$prtvtcfr== 2] <- 1
newdata$populist[newdata$prtvtfnl== 4|newdata$prtvtfnl== 3] <- 1
newdata$populist[newdata$prtvtbno== 8] <- 1
newdata$populist[newdata$prtvtdpl== 2 | newdata$prtvtdpl== 6] <- 1
newdata$populist[newdata$prtvtbse== 10] <- 1
newdata$populist[newdata$prtvtdes== 3|newdata$prtvtdes== 4 |newdata$prtvtdes== 6 |newdata$prtvtdes== 8| newdata$prtvtdes== 12] <- 1
newdata$populist[newdata$prtvtbgb== 7] <- 1
newdata$populist[is.na(newdata$populist)] <- 0


newdata$type[newdata$prtvtbat== 3] <- "populist-right"
newdata$type[newdata$prtvtcbe== 4|newdata$prtvtcbe== 7] <- "populist-right"
newdata$type[newdata$prtvtfch== 1] <- "populist-right"
newdata$type[newdata$prtvede1== 3] <- "populist-left"
newdata$type[newdata$prtvede1== 6] <- "populist-right"
newdata$type[newdata$prtvtdfi== 4] <- "populist-right"
newdata$type[newdata$prtvtbit== 4] <- "populist-left"
newdata$type[newdata$prtvtbit== 8|newdata$prtvtbit==  9] <- "populist-right"
newdata$type[newdata$prtvtcfr== 2] <- "populist-right"
newdata$type[newdata$prtvtfnl== 4] <- "populist-left"
newdata$type[newdata$prtvtfnl== 3] <- "populist-right"
newdata$type[newdata$prtvtbno== 8] <- "populist-right"
newdata$type[newdata$prtvtdpl== 2] <- "populist-right"
newdata$type[newdata$prtvtdpl== 6] <- "populist-right"
newdata$type[newdata$prtvtbse== 10] <- "populist-right"
newdata$type[newdata$prtvtdes== 3|newdata$prtvtdes== 4 |newdata$prtvtdes== 6 |newdata$prtvtdes== 8| newdata$prtvtdes== 12] <- "populist-left"
newdata$type[newdata$prtvtbgb== 7] <- "populist-right"
newdata$type[is.na(newdata$type)] <- "non-populist"

newdata$popl[newdata$prtvede1== 3] <- 1
newdata$popl[newdata$prtvtbit== 4] <- 1
newdata$popl[newdata$prtvtfnl== 4] <- 1
newdata$popl[newdata$prtvtdes== 3|newdata$prtvtdes== 4 |newdata$prtvtdes== 6 |newdata$prtvtdes== 8| newdata$prtvtdes== 12] <- 1
newdata$popl[is.na(newdata$popl)] <- 0

newdata$popr[newdata$prtvtbat== 3] <- 1
newdata$popr[newdata$prtvtcbe== 4|newdata$prtvtcbe== 7] <- 1
newdata$popr[newdata$prtvtfch== 1] <- 1
newdata$popr[newdata$prtvede1== 6] <- 1
newdata$popr[newdata$prtvtdfi== 4] <- 1
newdata$popr[newdata$prtvtbit== 8|newdata$prtvtbit==  9] <- 1
newdata$popr[newdata$prtvtcfr== 2] <- 1
newdata$popr[newdata$prtvtfnl== 3] <- 1
newdata$popr[newdata$prtvtbno== 8] <- 1
newdata$popr[newdata$prtvtdpl== 2] <- 1
newdata$popr[newdata$prtvtdpl== 6] <- 1
newdata$popr[newdata$prtvtbse== 10] <- 1
newdata$popr[newdata$prtvtbgb== 7] <- 1
newdata$popr[is.na(newdata$popr)] <- 0
summary(newdata$popr)

newdata$type2[newdata$prtvtbat== 3] <- "populist-right"
newdata$type2[newdata$prtvtcbe== 4|newdata$prtvtcbe== 7] <- "populist-right"
newdata$type2[newdata$prtvtfch== 1] <- "populist-right"
newdata$type2[newdata$prtvede1== 3] <- "populist-left"
newdata$type2[newdata$prtvede1== 6] <- "populist-right"
newdata$type2[newdata$prtvtdfi== 4] <- "populist-right"
newdata$type2[newdata$prtvtbit== 4] <- "populist-left"
newdata$type2[newdata$prtvtbit== 8|newdata$prtvtbit==  9] <- "populist-right"
newdata$type2[newdata$prtvtcfr== 2] <- "populist-right"
newdata$type2[newdata$prtvtfnl== 4] <- "populist-left"
newdata$type2[newdata$prtvtfnl== 3] <- "populist-right"
newdata$type2[newdata$prtvtbno== 8] <- "populist-right"
newdata$type2[newdata$prtvtdpl== 2] <- "populist-right"
newdata$type2[newdata$prtvtdpl== 6] <- "populist-right"
newdata$type2[newdata$prtvtbse== 10] <- "populist-right"
newdata$type2[newdata$prtvtdes== 3|newdata$prtvtdes== 4 |newdata$prtvtdes== 6 |newdata$prtvtdes== 8| newdata$prtvtdes== 12] <- "populist-left"
newdata$type2[newdata$prtvtbgb== 7] <- "populist-right"
newdata$type2[is.na(newdata$type)] <- "non-populist"

newdata$type2[newdata$prtvtbat== 2] <- "mainstream-right"
newdata$type2[newdata$prtvtcbe== 2|newdata$prtvtcbe== 3 |newdata$prtvtcbe== 9] <- "mainstream-right"
newdata$type2[newdata$prtvtfch== 4 |newdata$prtvtfch== 7] <- "mainstream-right"
newdata$type2[newdata$prtvede1== 1] <- "mainstream-right"
newdata$type2[newdata$prtvede1== 2] <- "mainstream-left"
newdata$type2[newdata$prtvtdfi== 3] <- "mainstream-right"
newdata$type2[newdata$prtvtbit== 1] <- "mainstream-left"
newdata$type2[newdata$prtvtbit== 5] <- "mainstream-right"
newdata$type2[newdata$prtvtcfr== 10] <- "mainstream-right"
newdata$type2[newdata$prtvtfnl== 2] <- "mainstream-left"
newdata$type2[newdata$prtvtfnl== 5|newdata$prtvtfnl== 7] <- "mainstream-right"
newdata$type2[newdata$prtvtbno== 7 | newdata$prtvtbno== 5] <- "mainstream-right"
newdata$type2[newdata$prtvtdpl== 1 | newdata$prtvtdpl== 4] <- "mainstream-right"
newdata$type2[newdata$prtvtbse== 3 | newdata$prtvtbse== 5] <- "mainstream-right"
newdata$type2[newdata$prtvtdes== 2] <- "mainstream-left"
newdata$type2[newdata$prtvtbgb== 1] <- "mainstream-right"
newdata$type2[is.na(newdata$type2)] <- "other-NA"

summary(newdata$populist)
newdata$type<-as.factor(newdata$type)
summary(newdata$type)
newdata$type2<-as.factor(newdata$type2)
summary(newdata$type2)

############separate variable for M5S (Italy)
newdata$m5s[newdata$prtvtbit== 4] <- 1
newdata$m5s[is.na(newdata$m5s)] <- 0

#################Select variables of interest and create a different dataset (ess)
myvars <- c("cntry","populist","ppltrst","pplhlp", "polintr", "trstprl", "lrscale", "stfeco","stfgov","gndr", "agea", "eisced", "ipcrtiv", "imprich",
            "ipeqopt", "ipshabt", "impsafe", "impdiff","ipfrule","ipudrst", "ipmodst","ipgdtim","impfree","iphlppl","ipsuces","ipstrgv", "ipadvnt",
            "ipbhprp", "iprspot", "iplylfr","impenv","imptrad","impfun", "type", "hinctnta","imueclt", "gincdif", "uempla", "uempli","type2","imwbcnt","imbgeco","imueclt", "sbstrec", "sbprvpv", "sbeqsoc", "sbbsntx", "sblazy", "sblwcoa", "trstep", "freehms", "hmsacld", "trstprt", "trstplt","pspwght","pweight","m5s","popl","popr")

ess <- newdata[myvars]

ess$income<-recode(ess$hinctnta,"77:100=NA")
ess$trustpeople<-recode(ess$ppltrst,"77:99=NA")
ess$interest<-recode(ess$polintr,"7:9=NA")
ess$trustpar<-recode(ess$trstprl,"77:99=NA")
ess$trustparties<-recode(ess$trstprt,"77:99=NA")
ess$trustpol<-recode(ess$trstplt,"77:99=NA")
ess$lr<-recode(ess$lrscale,"77:99=NA")
ess$economy<-recode(ess$stfeco,"77:99=NA")
ess$gov<-recode(ess$stfgov,"77:99=NA")
ess$gender<-recode(ess$gndr,"9=NA")
ess$age<-recode(ess$agea,"999=NA")
ess$education<-recode(ess$eisced,"55:99=NA;0=NA")
ess$highstudies<-recode(ess$eisced,"1:5=0;6:7=1;55:99=NA;0=NA")
ess$goodimmigration<-recode(ess$imueclt,"77:99=NA")
ess$tackleineq<-recode(ess$gincdif,"7:99=NA")
ess$unemployed[ess$uempla == 1] <- 1
ess$unemployed[ess$uempli == 1] <- 1
ess$unemployed[is.na(ess$unemployed)] <- 0

###########################################################recode basic human values 
ess$creative<-recode(ess$ipcrtiv,"7:9=NA")
ess$rich<-recode(ess$imprich,"7:9=NA")
ess$equality<-recode(ess$ipeqopt,"7:9=NA")
ess$admired<-recode(ess$ipshabt,"7:9=NA")
ess$safe<-recode(ess$impsafe,"7:9=NA")
ess$difference<-recode(ess$impdiff,"7:9=NA")
ess$rules<-recode(ess$ipfrule,"7:9=NA")
ess$diversity<-recode(ess$ipudrst,"7:9=NA")
ess$modest<-recode(ess$ipmodst,"7:9=NA")
ess$enjoy<-recode(ess$ipgdtim,"7:9=NA")
ess$free<-recode(ess$impfree,"7:9=NA")
ess$help<-recode(ess$iphlppl,"7:9=NA")
ess$sucess<-recode(ess$ipsuces,"7:9=NA")
ess$fist<-recode(ess$ipstrgv,"7:9=NA")
ess$adventure<-recode(ess$ipadvnt,"7:9=NA")
ess$behave<-recode(ess$ipbhprp,"7:9=NA")
ess$respect<-recode(ess$iprspot,"7:9=NA")
ess$loyal<-recode(ess$iplylfr,"7:9=NA")
ess$nature<-recode(ess$impenv,"7:9=NA")
ess$tradition<-recode(ess$imptrad,"7:9=NA")
ess$edonist<-recode(ess$impfun,"7:9=NA")
names(ess)

ess <- na.omit(ess)
ess[,c(76:96)]
ess[,c(76:96)]<- 6-ess[,c(76:96)]
summary(ess)


############################HOVs and centered HOVs
ess$security=rowMeans(ess[,c("safe", "fist")],)
ess$conformity=rowMeans(ess[,c("rules", "behave")],)
ess$tradition=rowMeans(ess[,c("modest", "tradition")],)
ess$benevolence=rowMeans(ess[,c("help", "loyal")],)
ess$universalism=rowMeans(ess[,c("equality", "diversity","nature")],)
ess$selfdirection=rowMeans(ess[,c("creative", "free")],)
ess$stimulation=rowMeans(ess[,c("difference", "adventure")],)
ess$hedonism=rowMeans(ess[,c("enjoy", "edonist")],)
ess$achievement=rowMeans(ess[,c("admired", "sucess")],)
ess$power=rowMeans(ess[,c("rich", "respect")],)

ess$mrat=rowMeans(ess[,c("safe", "fist", "rules", "behave", "modest", "tradition", "help", "loyal", "equality", "diversity", "nature", "creative", "free", "difference", "adventure", "enjoy", "edonist", "admired", "sucess", "rich", "respect")],)
ess$conformityc=ess$conformity-ess$mrat
ess$traditionc=ess$tradition-ess$mrat
ess$benevolencec=ess$benevolence-ess$mrat
ess$universalismc=ess$universalism-ess$mrat
ess$directionc=ess$selfdirection-ess$mrat
ess$stimulationc=ess$stimulation-ess$mrat
ess$hedonismc=ess$hedonism-ess$mrat
ess$achievementc=ess$achievement-ess$mrat
ess$powerc=ess$power-ess$mrat
ess$securityc=ess$security-ess$mrat

ess$openessc=rowMeans(ess[,c("hedonismc", "stimulationc", "directionc")],)
ess$enhancementc=rowMeans(ess[,c("achievementc", "powerc")],)
ess$conservationc=rowMeans(ess[,c("securityc", "conformityc", "traditionc")],)
ess$transcendencec=rowMeans(ess[,c("universalismc", "benevolencec")],)

########################anti-immigration and redistribution 
ess$immi1<-recode(ess$imbgeco ,"11:99999=NA")
ess$immi2<-recode(ess$imueclt ,"11:99999=NA")
ess$immi3<-recode(ess$imwbcnt ,"11:99999=NA")

ess$red1<-recode(ess$sbstrec ,"6:99999=NA")
ess$red2<-recode(ess$sbprvpv,"1=5;2=4;3=3;4=2;5=1;else=NA")
ess$red3<-recode(ess$sbeqsoc,"1=5;2=4;3=3;4=2;5=1;else=NA")
ess$red4<-recode(ess$sbbsntx ,"6:99999=NA")
ess$red5<-recode(ess$sblazy ,"6:99999=NA")

########################listwise NAs
myvars <- c("immi1", "immi2", "immi3","red1","red2","red3","red4","red5","universalismc","benevolencec","traditionc","securityc","conformityc","age","income","gender","highstudies","lr","cntry","type","type2","pspwght","pweight","m5s","populist","popl","popr")
ess <- ess[myvars]
ess<-na.omit(ess)
describe(ess)


#############factor anaysis immigration redistribution
myvars <- c("immi1", "immi2", "immi3","red1","red2","red3","red4","red5")
prueba <- ess[myvars]
prueba<-na.omit(prueba)

fit <- princomp(prueba, cor=TRUE)
summary(fit) # print variance accounted for
loadings(fit) # pc loadings
plot(fit,type="lines") # scree plot
fit$scores # the principal components
biplot(fit)

fit <- principal(prueba, nfactors=2, rotate="varimax")
fit

cfamodel <- ' immifac  =~ immi1 + immi2 + immi3
              redfac =~  red1 + red2 + red3 + red4 + red5'

fit <- cfa(cfamodel, data=prueba)
summary(fit, fit.measures=TRUE)
fitPredict <- as.data.frame(predict(fit))
ess <- cbind(ess, fitPredict)

#########################################weighting 
attach(ess)
ess$anweight<-pspwght*pweight
summary(ess)

####################################Table 3 model 1  (LPPs)
onlyleft <- ess[ which(ess$cntry =="DE" 
                       | ess$cntry == "IT"
                       | ess$cntry == "NL"
                       | ess$cntry == "ES"), ]

model1l<-glm(popl ~ universalismc + benevolencec + traditionc + securityc + conformityc + age + income + gender + highstudies + cntry, data = onlyleft, weights=anweight)
summary(model1l)
stargazer(model1l, type="text")
odd<-exp(coef(model1l))
stargazer(model1l, type="text", coef=list(odd), p.auto=FALSE)

####################################Table 3 model 2  (LPPs)
model2l<-glm(popl ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = onlyleft, weights=anweight)
summary(model2l)
stargazer(model2l, type="text")
odd<-exp(coef(model2l))
stargazer(model2l, type="text", coef=list(odd), p.auto=FALSE)

####################################Table 3 model 1  (RPPs)
onlyright <- ess[ which(ess$cntry != "ES"), ]
model1r<-glm(popr ~ universalismc + benevolencec + traditionc + securityc + conformityc + age + income + gender + highstudies + cntry, data = onlyright, weights=anweight)
summary(model1r)
stargazer(model1r, type="text")
odd<-exp(coef(model1r))
stargazer(model1r, type="text", coef=list(odd), p.auto=FALSE)

####################################Table 3 model 2  (RPPs)
model2r<-glm(popr ~ universalismc + benevolencec + traditionc + securityc + conformityc + age + income + gender + highstudies + redfac + immifac + lr + cntry, data = onlyright, weights=anweight)
summary(model2r)
stargazer(model2r, type="text")
odd<-exp(coef(model2r))
stargazer(model2r, type="text", coef=list(odd), p.auto=FALSE)


#############################Table 4: populist left VS mainstream left 
leftcontrol <- ess[ which(ess$type2 =="populist-left"| ess$type2 == "mainstream-left"),]
leftcontrol$mainstream[leftcontrol$type2 == "mainstream-left"] <- 0
leftcontrol$mainstream[leftcontrol$type2 == "populist-left"] <- 1

logistic<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + age + income + gender + highstudies + cntry, data = leftcontrol, weights=anweight)
logistic<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = leftcontrol, weights=anweight)
summary(logistic)
stargazer(logistic, type="text")
odd<-exp(coef(logistic))
stargazer(logistic, type="text", coef=list(odd), p.auto=FALSE)

###################### Figure OA.1 
plot(Effect("traditionc",logistic), confint=list(style="auto"), ylab="Probability of voting for Populist Left", xlab="Tradition", main=NULL)
plot(Effect("conformityc",logistic), confint=list(style="auto"), ylab="Probability of voting for Populist Left", xlab="Conformity", main=NULL)


####################################Table 4: populist right VS mainstream right 
rightcontrol <- ess[ which(ess$type2 =="populist-right"| ess$type2 == "mainstream-right"),]
rightcontrol$mainstream[rightcontrol$type2 == "mainstream-right"] <- 0
rightcontrol$mainstream[rightcontrol$type2 == "populist-right"] <- 1

logistic<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + age + income + gender + highstudies + cntry, data = rightcontrol, weights=anweight)
logistic<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = rightcontrol, weights=anweight)
summary(logistic)
stargazer(logistic, type="text")
odd<-exp(coef(logistic))
stargazer(logistic, type="text", coef=list(odd), p.auto=FALSE)


################################Figure OA.2
plot(Effect("universalismc",logistic), confint=list(style="auto"), ylab="Probability of voting for Populist Right", xlab="Universalism", main=NULL)
plot(Effect("traditionc",logistic), confint=list(style="auto"), ylab="Probability of voting for Populist Right", xlab="Tradition", main=NULL)
plot(Effect("securityc",logistic), confint=list(style="auto"), ylab="Probability of voting for Populist Right", xlab="Security", main=NULL)
plot(Effect("conformityc",logistic), confint=list(style="auto"), ylab="Probability of voting for Populist Right", xlab="Conformity", main=NULL)



####################################Last control: populist voting among ideologically moderate respondents (pg.14)
moderate <- rightcontrol[ which(rightcontrol$lr >=3 |rightcontrol$lr <=7),]
logistic3<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = moderate, weights = anweight)
summary(logistic3)

moderateb <- rightcontrol[ which(rightcontrol$lr <=7),]
logistic3b<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = moderateb, weights = anweight)
summary(logistic3b)

moderate2 <- leftcontrol[ which(leftcontrol$lr >=3 |leftcontrol$lr <=7),]
logistic5<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = moderate2, weights = anweight)
summary(logistic5)

moderate2b <- leftcontrol[ which(leftcontrol$lr >=3),]
logistic5b<-glm(mainstream ~ universalismc + benevolencec + traditionc + securityc + conformityc + immifac + redfac + age + income + gender + highstudies + lr + cntry, data = moderate2b, weights = anweight)
summary(logistic5b)